###############################################
########Loss-of-Strenght Gradient (OLS)########
###############################################

#Figure A.5

felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) +
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1","lag.lngcp",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.cost_sanction")])

omit<-NDs
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction+lag.lncapdist+lag.cost_sanction:lag.lncapdist+
                 lag.lnNTL+lag.lnpop_gpw_sum*lag.lncapdist+lag.lngcp*lag.lncapdist+lag.lnexcluded*lag.lncapdist
               +lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("LSG.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction", moderator = "lag.lncapdist",
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()
